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We have constructed a model for the kinetics of rupture of membranes under tension, apply- 
ing physical principles relevant to lipid bilayers held together by hydrophobic interactions. The 
membrane is characterized by the bulk compressibility (for expansion) K, the thickness 2ht of the 
hydrophobic part of the bilayer, the hydrophobicity a and a parameter 7 characterizing the tail 
rigidity of the lipids. The model is a lattice model which incorporates strain relaxation, and consid- 
ers the nucleation of pores at constant area, constant temperature, and constant particle number. 
The particle number is conserved by allowing multiple occupancy of the sites. An equilibrium "phase 
diagram" is constructed as a function of temperature and strain with the total pore surface and 
distribution as the order parameters. A first order rupture line is found with increasing tension, 
and a continuous increase in proto-pore concentration with rising temperature till instability. The 
model explains current results on saturated and unsaturated PC lipid bilayers and thicker artificial 
bilayers made of diblock copolymers. Pore size distributions are presented for various values of area 
expansion and temperature, and the fractal dimension of the pore edge is evaluated. 

PACS numbers: 87.16.Dg, 87.16. Ac, 87.14.Cc, 68.60.Dv, 05.10.Ln 



I. INTRODUCTION 

Fluid bilayer membranes made of lipids separate the 
contents of cells from their surroundings. The stabil- 
ity of those membranes under external perturbations is 
consequently of vital importance. In particular, if rup- 
tured the functions of the cell will be disabled. Natu- 
rally, a red cell bilayer membrane may rupture to free 
hemoglobin, protein responsible for oxygen transporta- 
tion in blood. This particular rupture, called blood cell 
hemolysis is achieved through thermal swelling of red 
cells |1j . Rupture or lysis can be achieved in a number of 
ways. Most experiments involved application of electric 
fields which produce compressive forces through the ca- 
pacitor effect. 01111 lllill. More recently a number 
of other ways have been used: mechanically, by suction 
through a pipette [j| , or extrusion through pores Q ; and 
by osmotic swelling of the cell \Hl IH E3 With the 
rapid progress in microscopic manipulation techniques, 
pore formation can also serve a positive purpose in drug 
delivery and gene therapy EE El, d The number of 
experimental studies of the mechanical properties of the 
cell is rapidly increasing and we only quoted a few exam- 
ples of recent work 01 • 

For the fluid bilayers, there is a significant degree of 
consensus on the order of magnitude of the quantities 
involved in membrane rupture. For phosphatidylcholine 
(PC) bilayers, area expansion seems to be only 2% to 4% 
before rupture. The corresponding external tension is of 
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the order of 1CT 3 to 1CT 2 N/m 0, S [3 E E3 . It also 
seems that in liquid membranes, rupture occurs via the 
nucleation of holes. Those observations are consistent 
with what is expected from structural considerations. 
Solid membranes as they go from brittle to ductile, expe- 
rience rupture kinetics which evolve from processes dom- 
inated by crack formation to hole formation 20j. More 
precisely, as the surface tension increases, the membrane 
will first be weakened by the apparition of small unstable 
pores with life-times as short as a few nanoseconds. At a 
critical tension, a large stable pore opens, relaxing almost 
entirely the membrane. This is rupture. Pore sizes can 
range from nano- to micro-meter in radius. In electropo- 
ration experiments, pores can last several microseconds 
0- 03- LH 0- EH- Using mechanical means, namely mi- 
cropipette extrusion on vesicles, pores can be kept open 
for several seconds before resealing [22, US • 

Several models of membrane rupture have been devel- 
oped over the years. Most of them derive from a model 
suggested by Litster a quarter of a century ago [24|- It 
explains the stability in terms of a surface energy and 
a pore edge energy. The model defines a critical pore 
size and an energy barrier for the creation of an irre- 
versibly growing pore. It has been extended and applied 
by other g roup s in particular to electric breakdown sit- 
uations [jj l25l |26| . Shillcock and Boal investigated the 
effect of temperature on membrane stability |27J . Tem- 
perature increases the entropy of the pore lowering its 
free energy. They show that even at zero tension, edge 
energy is required for stability. Netz and Schick have 
also developed a mean-field theory of the fluid bilayers 
as stacks of diblock copolymers [28| . 

In our work, we derive the energy of a finite size mem- 
brane under lateral expansion. We then consider pore 
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creation as a thermally activated process 29]. There is 
an important entropic element implied in a nucleation 
process. The model is for a bilayer whose relevant phys- 
ical quantities have their origin in the hydrophobicity of 
the lipid tails [3ll l3^ |. These quantities are the area 
compressibility (for expansion) K which results from the 
increased exposure of the hydrophobic tails as the mem- 
brane is stretched, the edge tension A, (or edge energy 
per unit length) which is the result of the exposure of 
these tails to water along the edge of the pore, and the 
rigidity 7 of the tails. 

In the model, the bilayer is characterized by an energy 
per site, essentially one molecule in size. It incorporates 
stress relaxation as a pore is created and grows. Pore- 
pore interactions are automatically taken into account. 
An equilibrium phase diagram is constructed determin- 
ing regions of no pores, non-critical pores (or protopores), 
and rupture. The phase diagram is in terms of the tem- 
perature and the area expansion of the membrane. The 
critical temperatures scale with the strength of the hy- 
drophobic interaction which to first order is linked to 
the length of the tails of the lipids. This hydrophobic 
interaction is obtained from parameters appropriate to 
pure phosphatidylcholine (PC) H El. Our model ap- 
plies to low strain rates since the model assumes that 
the membrane remains in mechanical equilibrium as it 
is stretched. We predict the rupture scenario for normal 
membranes, membranes made of unsaturated lipids, and 
very thick membranes. 

The following two sections develop the model. We first 
lay its physical foundations, and then present our Ising 
like model, which is solved using Monte Carlo simula- 
tions. We continue with the results and finish with a 
discussion of the possible extensions of the model. 
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where K is the compressibility (for expansion) and Aa = 
a — ao the excess area per molecule. For small expan- 
sion, the compressibility can be found experimentally 
from various techniques described previously. The ap- 
parent compressibility is the measure of the slope of the 
stress-strain curve of a membrane under a small tension, 
that is the slope of the curve 
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where r is the surface tension. However, since those 
membranes are fluid, thermal transverse undulations are 
usually present. For small stretches, a fraction of the ten- 
sion does not induce any strain, but rather brings back 
the membrane in its plane |33|. The real compressibil- 
ity is obtained by the slope of the stress-strain curve, 
but only in the regime where the undulations have been 
ironed out. For many lipids, this compressibility is about 
243 mN/m H3. 

When a single pore is inserted, some of this stress will 
be relieved. The relative change in the area per molecule 
becomes: 
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where a p is the area of the pore and Aa m is the total 
expansion of the membrane. But the edge of the new 
pore will now be exposed to water. As the tails are hy- 
drophobic, this will result in an increase of the energy 
proportional to the perimeter of the pore. For simplicity, 
we will assume the pore to be circular. In the presence of 
a pore, the energy of a stretched membrane then becomes 



II. MODEL FOR MEMBRANE RUPTURE 

Upon area increase, stress will build up in the mem- 
brane. As mentioned earlier, being a liquid the most 
likely stress release mechanism will be the formation of 
holes or pores. Let us imagine that a hole has appeared 
through thermal induced stress fluctuations. The ques- 
tion will be whether the gain in energy occasioned by 
the relaxation of the surfactants will counterbalance the 
energy loss through exposure of the tails of the molecule 
to the solvent. If this is the case, the hole will continue 
to grow into a large pore that will permit the system to 
relax almost entirely. For now, we assume the process 
to be purely planar. Whether fluctuations in the third 
dimension are important is an open question. At first 
thought they do not seem predominant. Assuming an 
elastic regime for small expansions of the membrane on a 
lattice of total relaxed area a m and molecular area ao, the 
energy cost associated with the stretching the membrane 
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The parameter A represents the edge tension (the effec- 
tive edge energy per unit length) . The first term in equa- 
tion J3J is the total edge energy, which represents the 
loss of energy on the perimeter of the circular hole. The 
second term, the surface energy, is the loss in energy re- 
sulting from the relaxation of the membrane induced by 
the hole, which augments the density of particles in the 
rest of the membrane. For the edge tension, it can be 
rewritten as 



A = 2h t <T, 



(5) 



where 2h t is the hydrophobic thickness of the bilayer, and 
therefore h t is the length of the lipid tails. And <r, the hy- 
drophobicity, is the hydrophobic energy of the lipid tails 
per unit area. With the usually quoted hydrophobicity 
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of 40 mN/ m J3llp3 . the edge tension A is of the order of 
1CP 7 mN [22,|32(. The experimentally measured values 
range from 1 to 4 x 10" 8 mN p^ . 

In the limit of large systems, the energy barrier to 
rupture is given by X 2 tt/ K(Aa m /a m )~ 1 and occurs for 
a pore area a p = (Xa m ^/Tr/K) 2 ^ 3 . When a p exceeds 
that critical value, the pore will grow till the mem- 
brane reaches the minimum energy configuration near 
a p = Aa rn . This argument predicts a barrier height of the 
order of W 4 k B T room (k B T room = 4.14 x lCT 18 mJ) for ex- 
perimentally observed stretches, such as Aa m /a m = 4% 
and a membrane of total area 1/im 2 , typical for PC mem- 
branes. Even if the cntropic components were very signif- 
icant as mentioned above, it is not sufficient to decrease 
significantly this large barrier. Moreover, according to 
these arguments, rupture should only happen at stretches 
around Aa m /a m = 50% ! 

In our opinion, the problem lies in an over-estimation 
of the line energy due to two factors. Firstly, when mem- 
branes expand, the tails get exposed to water, so ht is not 
a constant (especially not at 50%) . It decreases with in- 
creased stretching of the membrane. Secondly lipid bilay- 
ers are permeable and this permeability may be affected 
by the expansion of the membrane. So the hydrophobic- 
ity a used in the computation of the line energy will be 
smaller than the ideal value quoted above. We call a an 
apparent hydrophobicity and the a given above the pure 
hydrophobicity a pure . a will be an adjustable parameter. 



A. Compressibility, hydrophobicity and tail rigidity 

Under expansion, more water will enter the membrane, 
augmenting the free energy. The compressibility (upon 
expansion) K has therefore its origin in hydrophobic- 
ity. Moreover, the rigidity of the acyl tails controls the 
amount of water that penetrates the membrane at a given 
stretch. There is therefore a very close relationship be- 
tween surface energy, compressibility and tail rigidity. 
We will closely follow here some of the arguments of 
Wortis and Evans |32l (themselves based partly on Is- 
raelachvili's book |3lJ) for what pertains to the compress- 
ibility but adding the concept of the rigidity of the tails. 

To estimate K 7 we must understand the balance of 
forces that sets the area of the lipids a. The repulsive 
force between lipids is str ong at short distances, but falls 
off rapidly as a increases \?>'A . This repulsion is modelled 
by a potential of the form D/a, where D is a positive 
constant. On the other hand, the attractive part of the 
potential is a direct consequence of hydrophobicity. More 
water can enter in contact with the hydrocarbon chains 
as the lipids are separated from one another. In earlier 
work |3ll |32( a linear regime is assumed for small expan- 
sions of the cross section area of a lipid a, and the energy 
per molecule is written as a (a — <zo), where do is the 
equilibrium area. This assumes that the tails are totally 
flexible and that the area exposed to water is equal to the 
increase in area of the membrane (a — ag) . However, as 



illustrated in Fig. 1, the lipid tails more likely will only 
be able to partially close up the opening. We introduce 
a geometrical factor 7 > 1 defined so that j(a — eto) rep- 
resents the actual surface exposed to water per molecule 
under stretching, and hence (77(0 — ao) the increased en- 
ergy per molecule. The factor 7 is related to the rigidity 
of the tails, and therefore from now on, 7 will be referred 
as the tail rigidity. To have a geometric picture of 7, lets 
imagine that the stretching Aa is represented by a circle, 
and the effective exposed surface by a cone of surface S, 
as shown in figure 1. 7 is then the ratio S/ Aa. A rigid- 
ity of 1 represents tails that are completely flexible since 
the expanded lateral surface Aa is equal to the increased 
exposure of the tails to water. In contrast, a high rigid- 
ity represents very stiff chains, such that a small stretch 
can lead to a full exposure of the tails to water. With 
the geometry shown in Fig. 1, the length h e of the tails 
exposed to water can be written as : 



where r = yj a/n is the radius of the lipids. The non- 
exposed length h ne — h t — h e should be used in Eq. (J5J 
to calculate the line energy when a pore is nucleated at a 
certain stretch. Fixing the rigidity 7 for a certain type of 
lipid determines h e and h ne as a function of the effective 
stretch. 

If we combine both the attractive and the repulsive 
parts of the potential, we get a general expression for the 
energy per molecular site of the membrane: 

U(a) = 2<Tj(a~a ) + D(a- 1 ) + U Q . (7) 

The factor of 2 in the first term comes from the fact 
that we are dealing with a bilayer. The requirement 
dU/da\ ao = leads to the relationship D = 2a ja 2 ,. It 
must be emphasized that a here represents the interac- 
tion of the lipids when present in the membrane, and 
therefore refers to the apparent hydrophobicity, in oppo- 
sition to the pure hydrophobicity a pure discussed earlier 
in this section. By comparing the curvature of this po- 
tential around the minimum with the surface energy of 
Eq. (JTJ, we obtain 

K = 40-7. (8) 

Thus, there is a close relationship between the compress- 
ibility K, the hydrophobicity a and the rigidity 7. In 
our simulations, we use the experimentally known val- 
ues of compressibility K, and set the rigidity 7 to obtain 
rupture near 4% for the PC type bilayers. With the use 
of Eq. JHJ , we shall demonstrate in section IIVBI that 
the apparent hydrophobicity a is considerably smaller 
than the pure one cr pure , which ignores the water already 
present in the membrane. With a pure the line energy is 
over-estimated as we saw above. 
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III. THE CALCULATION 

What we are trying to simulate is a thermally activated 
nucleation process, it is therefore natural to consider a 
model amenable to a Monte Carlo simulation. For ease 
of computation we consider a lattice model similar to the 
well-known Ising model for binary mixtures, and consider 
the equilibrium phase diagram of the system as a func- 
tion of temperature and area expansion. This will reveal 
the expected scenarios of rupture as the membrane is 
expanded slowly, quasi-statically, so that the membrane 
remains in thermodynamic equilibrium at every step of 
the way. We will limit ourselves to two-dimensions, and 
to represent an isotropic liquid, we will use a hexagonal 
lattice (6 nearest neighbors). We begin by presenting 
the standard binary mixture model (SBMM) as applied 
to our problem and then describe the modifications and 
simplifications that we have made to it to incorporate 
the stress relaxation that occurs when a hole is created. 
In the SBMM, every site can be in either of two states, 
in our case "a" or "h" , representing respectively sites oc- 
cupied by a phospholipid or a hole. The site occupancy 
variable Si, is equal to 1 if the site is an "a" site and 
if it is a "h" site. The indices i and j run from 1 to 
N, where N is the total number of sites. Let J aQ (^ij), 
or its compact form J° a , be the interaction between two 



lipids separated by a distance fj 



IRi 



Rj |, where R, ; 



is the position vector for site i. In a similar way, J hh {rij) 
and J ah (rij) are defined. The microscopic Hamiltonian 
H consists of the sum of all interactions: 

H= \ E [JifsiSj + J^il-sMl-sj) 

z ij 

+ J^[*(l- Sj 0+*i(l-*)]] • (9) 

In a conventional canonical binary mixture, one would 
expand Eq. jSJ and eliminate the constant terms and 
the single sums since the number of each type of particle 
is conserved. In our binary mixture model, we keep the 
number N of particles (the a sites) and the total area 
constant. The number of holes (the h sites), however, 
can vary. Our model is neither the conventional canonical 
or grand canonical version of the binary mixture model. 
In the canonical ensemble, occupied (a) and empty (h) 
sites are interchanged to preserve the total number of 
particles. The conventional canonical model would be 
used if we were only interested in the phase separation of 
a system, which is not our case. In the grand canonical 
ensemble, the only dynamics that would be observed is 
the distribution of the holes. Neither model is suitable 
to incorporate the relaxation created by the appearance 
of a hole. 

We start with the SBMM in the canonical ensemble 
with an initial total number of particles equal to the num- 
ber of sites N. We now modify the model to impose the 
following constraints: the total area of the membrane 
stays constant as holes are created and the total number 
of lipids -a particles- stays constant. These lead to the 



following change in the dynamics. Occupied sites, when 
holes are created, now contain more than one particle, ac- 
tually N/(N — rih), where nt is the number of hole sites 
at a given time. The area of the occupied and hole sites 
stay constant. The nucleation of holes therefore relaxes 
the membrane and reduces the interparticle distance . 
This relaxation represents the surface energy and is the 
J°f term in Eq. ©. 

On the other hand, the interaction energy between 
holes and particles Jf- 1 represents the line energy of 
pores. The occupancy of sites being larger than one in 
such a case, there are more particles on the edge than the 
actual number of sites. As we will discuss in section llll Bl 
this will lead to a correction in the line energy. Finally 
the hole-hole interaction J-J- is set to zero. 



A. The surface energy 

The surface energy is not calculated using the first term 
in Eq. © but directly using the occupied area in the 
membrane. With nh hole sites, the interparticle distances 
have to be rescaled to the new values 
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where the m = | R, — Rj | are the initial intermolecular 
distances. As discussed earlier, the nucleation of holes 
changes the density of lipids in the membrane but con- 
serves their number. Since we assume uniform relaxation 
at each step of the simulation, the relative change in area 
per lipid, defined in Eq. Q, can be re- written as 
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where a p is now the total pore area, and not just the 
area of one pore. The surface energy can be calculated 
directly using the compressibility from the relation, 



E„ 
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It is clear from this last equation that if the relative 
number of holes n^/N completely relaxes the imposed 
stretch on the membrane Aa m /a m , the surface energy 
vanishes, as expected. 



B. The line energy 

As previously mentioned, the edge energy represents 
the hydrophobic contacts along the perimeter of a pore, 
measured by the hole- lipid interaction J^J 1 . The edge en- 
ergy therefore depends on the location of the holes in the 
lattice. The sum of those interactions cannot be reduced 
to a simple form like the surface energy. 
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As we are working on a hexagonal lattice, Jfp consists 
of the energy of exposing one sixth of the hydrophobic 
surface of a lipid to water, multiplied by two for the bi- 
layer. It is a first neighbor interaction. Jft can be writ- 
ten as 




if Tij is 1 



St. 



neighbor 
otherwise 



(13) 

where ao is the cross section of the lipids, and h ne = 
h t — h e the hydrophobic length of the lipid tails that are 
not exposed to water (ht and h e are defined in Eqs. (0 
and ®) respectively) . The hydrophobic energy of the 
exposed height h e is the origin of the surface energy. 

In the simulation the occupancy variables Sj remain 
equal to either or 1. So the edge energy as calculated 
by the program is 



T^modcl jah j ah 

^edge " ^ 1 
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where L ah is the number of lipid-hole interactions (be- 
tween sites). 

However, every time a hole site is created, the den- 
sity of the particles and hence the interparticle distance 
changes as expressed in Eq. <|10[l . in particular the near- 
est neighbor interparticle distance r . There are there- 
fore more particles on the edge of a pore than the actual 
number of sites. The edge energy per unit length can be 
written as J ah /f, where f is the scaled nearest neighbor 
interparticle distance. And hence the correct expression 
for the line energy is: 
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Using Eq. iJTUjl. Eq. (|ToT) can be written as 
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or in an expanded form, 
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C. Programming and Analysis 

To create and recover holes in the membrane up to 
an equilibrium point, we use the Metropolis algorithm in 
the Monte Carlo simulation. The mass of a hole is zero, 
but the mass of the system or number of particles is con- 
served through the renormalization of Eq. |JB}. Although 
the occupation of sites changes as the simulation evolves, 
detailed balance is satisfied because each configuration 
of pores has a uniquely defined energy. This energy does 
not depend upon the path followed to reach it. The tran- 
sition probability between a given initial and final config- 
uration will also be unique and reversible. If we choose 



an occupied site with the same probability than an empty 
site, that is, we select particles or holes according to their 
surface distribution; the acceptance rule will be given by 



acc(i -> /) = min (l,cxp - (js f - H^j /k B T ) , (18) 

in order to respect detailed balance. This acceptance 
rule is the usual one for a Monte Carlo simulation in the 
canonical ensemble. Finally, periodic boundary condi- 
tions are considered on the hexagonal lattice. 

The approach is highly efficient and large systems can 
be studied, up to 10 7 particles. However, systems of 10 4 
particles, which correspond to small vesicles, are nor- 
mally sufficient. Usually within 1000 samplings per site, 
the equilibrium state is easily reached. 



IV. RESULTS 

The types of membranes we study can be grouped in 
three main categories: saturated lipid bilayers (typical), 
unsaturated lipid bilayers, and long chain diblock copoly- 
mer bilayers. Membranes with unsaturated lipids may be 
found in biological systems. The long tail diblock copoly- 
mers, however, were synthesized by Bermudez et al., and 
then used to built vesicle shape polymersomes |3J|. In 
Table 1 are grouped the important parameters for our 
study. Those values are the total number of unsaturated 
bonds, the compressibility K, the membrane core thick- 
ness 2h t , the permeability, the critical tension for rupture 
t*, and the corresponding stretching Aa m */a m . Values 
for lipid membranes were taken from Olbrich et al. @ 
and Rawicz et al. 19]. The unsaturated lipid mem- 
branes differ from the typical membranes in the stiffness 
of their tails while the diblock copolymers simply have 
longer tails. 

In summary (details appear in the following subsec- 
tion) our strategy was to first use the above parameters 
to determine 7 the tail rigidity by looking at the effect it 
has on the critical rupture tension. K being known, every 
choice of 7 sets also the value of the apparent hydropho- 
bicity a through Eq. (jSJ) . Experimentally K seems inde- 
pendent of tail structure for most lipids [Tjj , but the rup- 
ture kinetics are not and this, we assume, is due mainly to 
changes in rigidity or 7. The compressibility of long tail 
copolymers is different from the usual 243 mN/m mainly 
because the he ad g roup is different, and not because of 
the tail length |34( • We will focus our attention mainly 
on the typical bilayers. 

Once 7 and hence a are set, we study the scenario 
of rupture through an "equilibrium" phase diagram that 
looks at the different qualitative changes in the distribu- 
tion of pores as a function of stretching for temperatures 
around kBT room . This corresponds to slow (quasistatic) 
stretching rates that keep the membrane in a state of 
equilibrium with respect to the applied stress. In the 
case of a higher rate of expansion, larger stretches and 
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TABLE I: Compressibility and other known physical param- 
eters of the different types of membranes studied. These pa- 
rameters are: the total number of cis-unsaturated bonds in 
the lipid tails ("uns."), the compressibility K , the thickness 
of the tails 2ht, the permeability P, the critical tension for 
rupture r*, and corresponding area stretching Aa m * /a m . 
Values for typical (i.e. with none or only one unsaturated 
bond) and unsaturated lipids are from Olbrich et al. || and 
Rawicz et al. |l9j|. whereas the values for long tail diblock 
copolymers are taken from Bermudez et al. |34|. 

uns. K 2h t P t* Aa m * /a r , 

membrane (mN/m) (nm) (fim/s) (mN/m) (%) 

Typical (average) 

0, 1 243±24 3,0±0,1 35±7 9±2 3,9±0,9 1 

Unsaturated lipids 

SLPC:0/2 2 243±24 3,0±0,1 49±6 4,9±1,6 2,0±0,5 1 
DLPC:2/2 4 243T24 3,0±0,1 91±24 5,1±1,0 2,1±0,4 1 
DLPC:3/3 6 243±24 3,0±0,1 146±24 3,1±1,0 1,2±0,4 1 
Long tail copolymers 

OE7 102±10 8±1 21±3 1 21±2 

OB9 102±10 11±1 28±4 1 28±3 

OB18 102±10 21±1 — 40±5 1 40±4 
Abbreviations: SLPC:0/2 = 
l-stearoyl-2-linoleoyl C isat9,i2-phosphatidylcholine; 
DLPC:2/2 = dilinoleoyli, ot ( lci3at g i i2-phosphatidylcholine; 
DLPC:3/3 = dilinoleoyl i , othci3at 9 i i2,i5-phosphatidylcholine; 
OE7 = ethylenoxide4o-ethylethylene37; 
OB9= ethylenoxideso-butadieness; 
OB18= ethylenoxide8o-butadienei25 . 
^Values obtained using Eq. ||5J, and strains or stresses to 
rupture from quoted articles. 



therefore tensions will be needed to break the membrane 
[29l l30| . Because we are considering thermally activated 
processes, higher temperatures in this model are equiv- 
alent to weaker interactions. To study the reversibility 
of membrane lysis, hysteresis curves were constructed to 
compare the number of holes and their distribution as 
we first stretch the membrane to rupture, and then let it 
relax to zero tension. This study revealed that the rup- 
ture transition is first order at room temperature, with a 
strong hysteresis. The fractal dimensions of stable pores 
were computed at different temperature as a measure of 
their entropy and shape. As a prelude we look at system 
size dependence. 



A. System size dependence 

To simplify the discussion, we assume that the correct 
tail rigidity has been obtained. In reality, the numerical 
value of the tail rigidity depends somewhat on system 
size. 

Bilayer membranes can form vesicles of a variety of 
shapes. In nature, they are normally present as spherical 
vesicles of radii between 1 to 10/im [11-13]. For a given 
relative stretch Aa m /a m the surface energy per particle 
does not depend on the size of the lattice. The line energy 
scales as the length of pore edge or the sum of the square 



root of the surfaces of the individual pores. System size 
may affect the kinetics because of the changes in pore 
size distribution and the fact that pore size distribution 
does not scale with system size. 

As shown in figure 2, stretching to rupture almost sta- 
bilizes to about 3.9% for systems larger than 10 4 par- 
ticles. For smaller systems, we observe a significant in- 
crease. As expected, this plateau is mainly due to the 
increase in the number of pores at rupture, which raises 
the line energy (see figure 3). Furthermore, entropy in- 
creases for larger systems, which favors earlier rupture. 

.These plots are averages with a standard deviation of 
10 runs. It is an interesting and open question to know 
if multiple pore rupture is observed in biological mem- 
branes. It must be noted that the number of sites in our 
simulations corresponds to half the number of particles, 
since we are dealing with bilayer membranes. For what 
follows, we shall work on a lattice of 30301 sites (60602 
lipids). This corresponds to a total membrane area of 
18/im 2 (ao ~ 0.6nm 2 for phospholipids). For a spheri- 
cal vesicle, it corresponds to a radius of 1, 5/im, a small 

"natural biological cell. 



B. The rigidity and hydrophobicity of lipid tails 

Depending on the rigidity of the tails of a lipid, more or 
less water may penetrate the membrane under stretching. 
For unsaturated lipids, one would expect higher rigidi- 
ties, since the tails are stiffer. Furthermore, if the tails 
are rigid, the membrane is also more permeable |19| . In 
this case, the apparent hydrophobicity of the membrane 
is lower, since the difference in energy for the lipid inside 
and outside the membrane is smaller. For small varia- 
tions of these parameters, we can see from Eq. 7 that 
a saturated and an unsaturated membrane made out of 
the same lipid group should have comparable compress- 
ibilities. This was verified experimentally by Rawicz et 
al. 0| in their investigation of the effects of lipid un- 
saturation on membrane elasticity. This group has found 
that unsaturation has no effect on compressibility. Never- 
theless, unsaturation is still an important issue for mem- 
brane rupture, as unsaturated lipids are less flexible than 
saturated ones. Olbrich et al. § have discovered that 
rupture tension was only affected for bilayers made with 
lipids of two or more alternating cis-double bonds (C=C- 
C=C) in one or both chains. The three degrees of un- 
saturation that we study are as shown in Table 1: lipids 
with one tail having two alternating cis-double bonds (1 
pair), both tails having two alternating cis-double bonds 
(2 pairs) and both tails with three alternating cis-double 
bonds (4 pairs). All these unsaturated lipids also have 
18 carbon atoms in each tail. For the long polymers, the 
line energy is directly proportional to the length of the 
lipid tails. It is therefore natural to expect higher ten- 
sions at rupture for longer chains. Membrane rupture of 
diblock copolymers membranes will also be studied. 

To obtain 7 and therefore a, of the different mem- 



7 



TABLE II: Tail rigidity and apparent hydrophobicity of the 
different lipids studied. The estimated value for pure hy- 



drophobicity is also given 


as a 


_ _ r„ _ _ __ „_ 

rcicrciicc. 






Membrane 


Name 


Uns. 




(J 












(mN/ 


m) (mN/m) 


Typical 


(average) 


0, 1 


8.0±1.0 


7.6 


40 


Unsaturated 


SLPC:0/2 


2 


13.0±3.0 


4.7 


w 40 


lipids 


DLPC:2/2 


4 


13.0±3.0 


4.7 


w 40 




DLPC:3/3 


6 


17.5±0.5 


3.8 


w 40 


Long tails 


0E7 





7.3±1.0 


3.4 




copolymers 


OB9 





8.6±1.0 
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13.0±1.0 


2.0 





branes, we fit the tail rigidity to have the correspond- 
ing stretching to rupture. For typical lipid bilayer mem- 
branes, we need a rigidity of 7 = 8 ± 1 to have rupture 
at 3.9% (see figure 4). From Eq. (7), this corresponds to 
a hydrophobicity of 7.6. 

The plateau near 2% stretching in figure 4 for large 7's 
may seem abnormal at first sight. Even more, just be- 
fore zero tension rupture, there is even a slight increase in 
the stretching to rupture. This effect is related to the cri- 
terion that defines rupture. Zero tension rupture means 
that the lipids are so stiff and the membrane so permeable 
that stability is not even assured. Just before this critical 
rigidity, the membrane is still not very stable, since the 
line energy per lipid rapidly decreases with stretching. 
This weak interaction between lipids raises the entropy 
of the system. For high rigidities, this gain in entropy will 
make the rupture scenario more comparable to a melting 
phenomenon than a pore creation relaxation. Since we 
have lost the ideal scenario of rupture which favors relax- 
ation as opposed to the creation of line energy, rupture 
in terms of a few stable pores is more difficult to achieve. 
The definition of rupture will be discussed in more detail 
in the next subsection. 

Nevertheless, this plateau seems consistent with ex- 
perimental observations. Unsaturated lipids with 2 and 
4 total unsaturated bonds both have a stretching to rup- 
ture near 2% (see figure 4). This therefore means that 
even though they have the same stretching to rupture 
does not mean they have the same tail rigidity. In fact, 
they should not. 

For typical membranes, the apparent hydrophobicity 
a is about 5 times smaller than the approximated value 
for (jp Ure (see table 2). This leads to a line energy of 
2.4 x 10 -8 mN, well with the experimentally measured 
range of 1 to 4 x 10 -8 mN. For highly unsaturated lipids, 
this ratio can be as high as 10 if we suppose that the pure 
hydrophobicity has not changed significantly. It would 
be wrong, however, to say that 80% of the membrane 
is flooded. This apparent hydrophobicity may also take 
into account other repulsive energies, such as the head- 
head repulsion. Moreover, a pure was approximated from 
geometrical simplifications, and might be smaller than 40 



mN/m. Nevertheless, this noticeable difference between 
the pure and the apparent hydrophobicity indicates that 
membrane permeability is not a negligible effect. 

Graphics similar to figure 4 were used to obtain the tail 
rigidity of the diblock copolymers studied by Bcrmudcz 
et al. ;341] . The tail rigidity of OE7 and OB8 is compa- 
rable to the rigidity of a typical lipid tail. However, for 
OB18, which has a much longer tail (see Table I) 7 is 
considerably higher. It may be that for very long poly- 
mers tail entanglements make the tails appear more rigid 

mm. 

In the following subsections, we focus on typical bio- 
logical membranes. So 7 is set equal to 8. 



C. Phase diagram 

A phase diagram was constructed to study the behav- 
ior of the membrane under variations of the temperature 
and the surface area (i.e. under stretching). We show 
the phase diagram in Fig. 5 for a "typical" membrane 
(see table I). The most striking feature is a "first order 
transition" line to a ruptured state. When at a given tem- 
perature, the membrane is stretched, one observes in the 
neighborhood of a specific stretch an abrupt increase in 
the relaxation which we call rupture. This occurs through 
the formation of a few massive stable pores. The line 
shown in Fig. 5 corresponds to the point of abrupt re- 
laxation as the tension is increased. Below this line the 
nature of the membrane evolves gradually with tempera- 
ture T. At low T, including around T room , there is what 
we call a stable state. The relaxation in this regime is 
small, with occasionally a few protopores (see Section 
II V E|) . As the temperature is raised the number and the 
size of the protopores keep increasing. The pores become 
numerous, but are still short lived, and of size no larger 
than 10 particles or so. 

At temperatures below and near T room , rupture is very 
sharp. We go directly from a state of high tension where 
less than 20% of the membrane is relaxed to an almost 
fully relaxed system (80% and more) . At higher temper- 
atures, there are a larger number of pores present in the 
membrane at low stretches, but they do not aggregate 
till the tension reaches the rupture tension. 

The rupture line first decreases with temperature till 
1.6T room then rises again. Entropy appears to favor a 
high density of dispersed pores. It takes then an in- 
creased tension to force the aggregation of the pores and 
trigger the full relaxation of the membrane. This leads 
to the rise in the rupture tension. Finally the rupture 
line ends near T = 2.8T room , but it is no longer a first 
order transition line. At that temperature the membrane 
loses stability: stable pores are created without any ten- 
sion applied. The reduced line energy A/fc^T falls below 
the critical value required to assure membrane stability. 
Shillcock and Boal [27] vary A in their tethered beads 
model. They find a critical value of A*6/fcsT = 1.66 at 
constant zero-tension, where b is the average vertex sepa- 
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ration in their model (At T room , this means a critical line 
energy of about 0.9 x 10~ n J/m). We work at constant 
area and, in our model, we calculate hole-lipid interac- 
tions and the total line energy is given by J ah L ah (see 
Eq. JUJl) , where J?- is the hole-lipid interaction defined 
in Eq. (|13|) . The effective line tension depends on pore 
distribution. With our thermally activated model the be- 
havior observed near 2.8T room can be viewed as behavior 
near T room with a line energy reduced by the factor 2.8. 
Using the distribution of pores in Fig. 0:, one can find 
an effective line energy per site of approximately l.U ah 
(small pores dominate the distribution) j3|| , or a critical 
line energy per unit length of about 0.94 x 10~ n J/m. 
Interestingly our lattice model and Shillcock and Boal's 
bond flipping model have consistent views on stability. 
Mechanically the membrane may fail earlier so we label 
the region around T — 2.8T room as unstable. 

At room temperature, the presence of a regime offer- 
ing 2% to 20% unstable pore surface, or small non-stable 
pores, called protopores, have experimentally been ob- 
served at room temperature p|. In the next three sub- 
sections we look at different properties of the membrane 
under tension that give a clearer picture of the rupture 
kinetics. 



D. Hysteresis 

To show that indeed rupture is "first-order", samples 
were put through a full cycle of expansion and compres- 
sion. This is done again in a quasi-static process, i.e. at 
each step the system is allowed to relax. We follow, as the 
strain was varied, the variations in the stress (Fig. 6a), 
the relative pore area (a measure of the relaxation) (Fig. 
6b), and the total number of pores (a measure of pore 
interaction and coalescence) (Fig. 6c). 

Hysteresis effects are found to be very strong at low 
temperature. Rupture and healing occur through nuclc- 
ation of one pore. At high temperature, hysteresis is 
reduced and there are many more pores present, below 
and above rupture. At T roomi the case illustrated in Figs. 
6, behavior is similar to the high temperature situation 
below rupture, but comparable to the low temperature 
case above rupture. A considerable number of pores are 
present in the membrane at first, but after rupture, only 
a small number of large non-fluctuating pores remain. 
This is clearly observable in Fig. 6c as the rapid rise in 
the number of pores which precedes the collapse at rup- 
ture. The tension at collapse defines the T room point on 
the rupture line. 



E. Pore size distribution 

To further illustrate the difference in pore distribution 
below and above rupture, we compare the distribution of 
pore sizes for membranes with 2% and 6% strain. We do 
this at the three temperatures T room , 2T room , and 3T room 



to illustrate the change in structure of the membrane 
with increasing temperature (see Figs. 13(a), (b), and (c)). 
The distributions are typical snap-shots of the membrane 
as the Monte Carlo simulation evolves. For the stable 
membranes shown, we see in figures 7a and 7b, the sharp 
drop in the number of small pores above rupture. A large 
pore creates less line energy than a protopore for the 
same gain in surface energy, so is energetically favored. 
For unstable membranes, the appearance of larger pores 
in the system does not reduce the number of smaller ones 
(figure 7c). The small pores, present at zero tension at 
3T room are stable due to the contribution of entropy to 
the free energy. This illustrates from another perspective 
the instability of the membrane at 3T room . 



F. Pore shape 

At different temperatures for a stretching of 3% (just 
before rupture), we insert into the membrane a hole at 
the origin. After the pore is fully opened and stable, 
we compute its fractal dimension df from the scaling re- 
lationship between the pore area A and its perimeter I 
(A = A(0)l d f where A$ is a proportionality factor). 

From the results plotted in figure 8, we clearly see that 
df decreases as the temperature is raised. This is caused 
by increased irregularities along the edge of the pore. 
At T = 2.8 T room , where rupture occurs spontaneously 
at zero tension, one would expect the scaling relationship 
applicable to a self-avoiding ring [27| ■ In such geometries, 
the area of the pore should scale as A — AqI 3 / 2 . For this 
case, our model gives df — 1, 54 ±0, 10 in agreement with 
this argument. 

V. CONCLUSION 

Bilayer membranes are an intrinsic part of all living 
species. With the rapid development of biotechnology, 
membranes can now be artificially created, and there- 
after even modified to form vesicles of controlled size and 
properties. Rupture kinetics play an integral part in the 
modifications of those vesicles and in their applications 
such as in drug delivery. 

In spite of their complex molecular structure, lipid 
bilayer membranes offer an ideal model system, as hy- 
drophobicity is responsible for both the edge line energy 
and the bulk compressibility (for extension). We have 
developed the first model for the rupture of a membrane 
held together by hydrophobic forces, which includes the 
nucleation and growth of pores. Our minimal model can 
be used to reproduce current results on saturated and 
unsaturated PC lipid bilayers and thicker artificial bilay- 
ers made of diblock copolymers. We have introduced a 
new quantity called the rigidity of the tails 7 which per- 
mits the study the effect of saturation of the lipid tails 
on the structural properties of the bilayer. Interestingly 
an increase in rigidity produces a decrease in apparent 
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hydrophobicity since K is nearly constant for a given 
class of lipids. This may mean an increase in permeabil- 
ity. The model requires little computer time allowing the 
handling of real size vesicles. Entropic and nucleation ef- 
fects are included naturally. 

The structural integrity of a lipid bilayer is consid- 
erably affected by its composition; biological membrane 
structure can be fairly complex. Different types of lipids, 
cholesterols, and proteins play an important role in deter- 
mining membrane stability. The model has the potential 
of handling these inclusions. 
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FIG. 1: Schematic representation of the total height of the 
lipid tails exposed to water at a given expansion. The rigidity 
7 is the ratio of surface of the sides of the cones to the base. 
The total surface exposed to water is expected to be larger 
than the surface expansion. 
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FIG. 2: The dependence of the degree of stretching on the 
size of the membrane for rupture at room temperature. The 
number of sites is half the number of particles since we are 
working with a bilayer. This value stabilizes at area expan- 
sions near 4% for systems larger than 10 4 sites. The dashed 
line indicates the size of the network used in the rest of the 
simulations. 
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FIG. 3: The average number of stable pores in the system 
when it ruptures at room temperature. The relation is linear 
(exponential on this log scale). The dashed line indicates the 
size of the network used in the rest of the simulations. 
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FIG. 4: The degree of stretching for rupture at room tem- 
perature for different tail rigidities for membranes with com- 
pressibility of 243 mN/m and membrane core thickness of 3,0 
nm (see table 1). The dashed lines are the known area ex- 
pansions at rupture for membranes with different degrees of 
unsaturation. 
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FIG. 5: Phase diagram predicted by the model for a "typical" 
bilayer (se Table I) showing the degree of relaxation in the 
membrane as a function of stretching and temperature. The 
full square dots indicate the rupture first-order transition line. 
The dashed lines are the curves of constant relative relaxation. 





FIG. 6: Hysteresis of different quantities upon stretching at 
T = T room : (a) the tension r; (b) the relative area occupied by 
pores which is equal to the number of pore sites divided by the 
total number of sites; (c) the total number of pores. Room 
temperature systems behave similarly to high temperature 
systems at low tensions by having a considerable quantity of 
protopores present in the membrane. The rupture scenario is 
however similar to that of lower temperature systems. 
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FIG. 7: Distribution of pores at three temperatures (a) 
T = Troom, (b) T — 2T room , and (c) T = 3T room , below 
and above the rupture point in a lattice of 30301 sites (each 
site « 0.6 nm 2 ). With increasing temperature, the number 
of small pores increases and a larger number remain in the 
system after rupture. At T = 3T room , above the instability 
point of 2.8Troom, there is little structural difference between 
a 2% and 6% stretched membrane except for the larger size 
of the pores. 
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FIG. 8: The Fractal dimension of a pore in a stretched mem- 
brane. This quantity is a measure of the regularity of the 
shape of the pore and its edge. As the temperature is in- 
creased, the pore becomes irregular and spreads over a larger 
membrane surface. 



